/*
 * interaction_4HSiC.cpp
 *
 *  Created on: 2011-11-9
 *      Author: nzhao
 */


#include "special/interaction_4HSiC.h"

double  * NE_INT_MAT1, * NE_INT_MAT2;  //interactions between electron spin and bath spins
double	* MEDIATED_INT1,* MEDIATED_INT2;
double	* MAX_MEDIATED_INT;

double 	NZEEMAN;
double	* HF_TENSOR;

double 	EIGEN_VALUES[3];
double	S1VALUE[3], S2VALUE[3];
double	bframeX[3], bframeY[3], bframeZ[3];


int			count_core;
double		** RMatrix;
CLD			** SMatrix;
CLD 		** S_ab;
CLD 		** S_ab_Bframe;

CLD SX[9] ={CLD(0.,0.),			CLD(1/SQRT2, 0.0), 	CLD(0.0, 0.0),
			CLD(1/SQRT2,0.),	CLD(0.0, 0.0), 		CLD(1/SQRT2, 0.0),
			CLD(0.,0.),			CLD(1/SQRT2, 0.0), 	CLD(0.0, 0.0)};

CLD SY[9] ={CLD(0.,0.),			CLD(0.0, -1/SQRT2), 	CLD(0.0, 0.0),
			CLD(0.0, 1/SQRT2),	CLD(0.0, 0.0), 		CLD(0.0, -1/SQRT2),
			CLD(0.,0.),			CLD(0.0, 1/SQRT2), 	CLD(0.0, 0.0)};

CLD SZ[9] ={CLD(1.,0.),	CLD(0.0, 0.0), CLD(0.0, 0.0),
			CLD(0.,0.),	CLD(0.0, 0.0), CLD(0.0, 0.0),
			CLD(0.,0.),	CLD(0.0, 0.0), CLD(-1.0, 0.0)};

void	set_b_frame();
void	uni_transformation(int dim, CLD * mat, CLD * R_mat, CLD * res_mat);
void 	Crystal2Nature(CLD * coord_c, CLD * coord_n);
void 	Nature2Bframe(CLD * coord_n, CLD * coord_b);

double 	vec_dot(int dim, double * vec1, double * vec2);
void 	vec3_cross(double * vec1, double * vec2, double * vec3, bool is_norm);
void 	vec_dot_mat(CLD * S_ab, double * hftensor_i, CLD * sab_A_vec);

void calc_Interactions()
{
	int i, j, k;

	set_b_frame();

	CentralEigen( THETA, PHI);

	/////////////////////////////////////////////////////////////////////////
	// interactions for Hami1
	NE_INT_MAT1      = new double [NBATH * 3];
	NE_INT_MAT2      = new double [NBATH * 3];

	for (i = 0; i < NBATH * 3; i++)
	{
		NE_INT_MAT1[i] = 0.0;
		NE_INT_MAT2[i] = 0.0;
	}

	count_core = 0;
	for (i = 0; i < NBATH; i++)
	{
		calc_hf_1(i, NE_INT_MAT1);
	}

	count_core = 0;
	for (i = 0; i < NBATH; i++)
	{
		calc_hf_2(i, NE_INT_MAT2);
	}

	//////////////////////////////////////////////////////////////////////////
	//interactions for mediated nuclear spin interactions
	if(IsMediated)
	{
		HF_TENSOR 		= new double [NBATH * 9];
		MEDIATED_INT1	= new double [NBATH * NBATH * 9];
		MEDIATED_INT2	= new double [NBATH * NBATH * 9];
		MAX_MEDIATED_INT = new double [NBATH * NBATH];

		for(i = 0; i < NBATH; i ++)
			calc_hf_tensor(i);

		int max_count = 0;
		for(i = 0; i < NBATH; i ++)
		{
			for(j = 0; j < NBATH; j++)
			{
				calc_mediated_int(i, j, STATE_1, MEDIATED_INT1 + (i * NBATH + j) * 9 );
				calc_mediated_int(i, j, STATE_2, MEDIATED_INT2 + (i * NBATH + j) * 9 );

				MAX_MEDIATED_INT[i * NBATH + j] = 0.0;
				for(k = 0; k < 9; k++)
				{
					if( fabs(MEDIATED_INT1[ (i * NBATH + j)*9 +k ]) >  MAX_MEDIATED_INT[i * NBATH + j])
						MAX_MEDIATED_INT[i * NBATH + j] = fabs(MEDIATED_INT1[ (i * NBATH + j)*9 +k ]);
					if( fabs(MEDIATED_INT2[ (i * NBATH + j)*9 +k ]) >  MAX_MEDIATED_INT[i * NBATH + j])
						MAX_MEDIATED_INT[i * NBATH + j] = fabs(MEDIATED_INT2[ (i * NBATH + j)*9 +k ]);
				}
				if(MAX_MEDIATED_INT[i * NBATH + j] > MEDIATED_INT_THRESHOLDVALUE)
					max_count ++;

			}
		}
		cout << max_count / 2 << "pairs exceed the interaction threshold." << endl;
	}
}

void CentralEigen(double t, double p)
{
	int i, j, k;


////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
// Central spin Hamiltonian
	CLD central_hami[9], evec[9];
	double eval[3];

	double bx, by, bz;
	bz = - EXT_MAG_FIELD * GAMMA_E / 2. / PI * cos(t) / 1.E9;
	bx = - EXT_MAG_FIELD * GAMMA_E / 2. / PI * sin(t) * cos(p) / 1.E9;
	by = - EXT_MAG_FIELD * GAMMA_E / 2. / PI * sin(t) * sin(p) / 1.E9;

	for(i = 0; i < 9; i++)
	{
		central_hami[i] = bx * SX[i] + by * SY[i] + bz * SZ[i];
	}
	central_hami[0] += ZFS_4H_SiC;
	central_hami[8] += ZFS_4H_SiC;

	MKL_hermite_mat_diag(3, central_hami, eval, evec);

	for(i = 0; i < 3; i++)
		EIGEN_VALUES[i] = eval[i] * 2.0 * PI * 1.E6;// GHz to k rad / s

	cout << "eigen values (in GHz):" << endl;
	cout << eval[0] << "\t" << eval[1] << "\t" << eval[2] << endl;

////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
// S_alpha_beta matrix
	CLD SXMatrix1[9], SYMatrix1[9], SZMatrix1[9];

	uni_transformation(3, SX, evec, SXMatrix1);
	uni_transformation(3, SY, evec, SYMatrix1);
	uni_transformation(3, SZ, evec, SZMatrix1);

	CLD vec[3];

	S_ab 		= new CLD * [3 * 3];
	S_ab_Bframe = new CLD * [3 * 3];
	RMatrix = new double * [3 * 3];
	SMatrix = new CLD * [3 * 3];

	for(i = 0; i < 9; i ++)
	{
		RMatrix[i] = new double [3 * 3];
		SMatrix[i] = new CLD [3 * 3];
		S_ab[i] = new CLD[3];
		S_ab_Bframe[i] = new CLD[3];
	}

	int index_ab, i1, j1;
	for(i = 0; i < 3; i++)
	{
		for(j = 0; j < 3; j ++)
		{
			index_ab = i * 3 + j;
			vec[0] = SXMatrix1[index_ab]; vec[1] = SYMatrix1[index_ab]; vec[2] = SZMatrix1[index_ab];

			Crystal2Nature(vec, S_ab[index_ab]);
			Nature2Bframe(S_ab[index_ab], S_ab_Bframe[index_ab]);

			//Mean spin value
			if(i == STATE_1 && j == STATE_1)
			{
				S1VALUE[0] = S_ab[index_ab][0].real();
				S1VALUE[1] = S_ab[index_ab][1].real();
				S1VALUE[2] = S_ab[index_ab][2].real();
			}
			if(i == STATE_2 && j == STATE_2)
			{
				S2VALUE[0] = S_ab[index_ab][0].real();
				S2VALUE[1] = S_ab[index_ab][1].real();
				S2VALUE[2] = S_ab[index_ab][2].real();
			}

			//outer product to generate a tensor
			for(i1 = 0; i1 < 3; i1 ++)
			{
				for(j1 = 0; j1 < 3; j1 ++)
				{
					SMatrix[index_ab][i1 * 3 + j1] = S_ab_Bframe[index_ab][i1] * conj(S_ab_Bframe[index_ab][j1]);
					RMatrix[index_ab][i1 * 3 + j1] = 2.0 * SMatrix[index_ab][i1 * 3 + j1].real();
				}
			}
		}
	}

	cout << "S1 value:" << STATE_1 << "\t" << "S2 value:" << STATE_2 << endl;
	cout << setw(12) << setiosflags(ios::scientific) << setprecision(10);
	cout << S1VALUE[0] <<"\t" << S2VALUE[0] << endl;
	cout << S1VALUE[1] <<"\t" << S2VALUE[1] << endl;
	cout << S1VALUE[2] <<"\t" << S2VALUE[2] << endl;
	cout << endl;
}

///////////////////////////////////////////////////////////////////////////////////////////////
//intrinsic dipolar interactions
sINT_TENSOR1 calc_dipole_int(int idx1, int idx2)
{
	int i, idxij;
	double * d = new double [3];
	double * n = new double [3];
	double len, len3, nx, ny, nz;

        double int_const_n1n2;
        int_const_n1n2 = INT_CONST_NN * RELATIVE_GAMMA_N[idx1] * RELATIVE_GAMMA_N[idx2];

	sINT_TENSOR1 res;

	idxij = 6 * (idx1 * (idx1 + 1) / 2 + idx2);

	if(idx1 != idx2)
	{
		len = 0.;
		for (i = 0; i < 3; i++)
		{
			d[i] = COORDINATES[3 * idx1 + i] - COORDINATES[3 * idx2 + i] ;
			len += d[i]*d[i];
		}
		len = sqrt(len);
		len3 = len * len * len;

		for (i = 0; i < 3; i++)	{	n[i] = d[i] / len; }

		nx = vec_dot(3, n, bframeX);
		ny = vec_dot(3, n, bframeY);
		nz = vec_dot(3, n, bframeZ);


		res.xx = int_const_n1n2 * (1.- 3.* nx * nx) / len3;
		res.yy = int_const_n1n2 * (1.- 3.* ny * ny) / len3;
		res.zz = int_const_n1n2 * (1.- 3.* nz * nz) / len3;
		res.xy = int_const_n1n2 * (  - 3.* nx * ny) / len3;
		res.yz = int_const_n1n2 * (  - 3.* ny * nz) / len3;
		res.zx = int_const_n1n2 * (  - 3.* nz * nx) / len3;
		res.yx = res.xy;
		res.zy = res.yz;
		res.xz = res.zx;
	}
	else
	{
		res.xx = 0.0;
		res.yy = 0.0;
		res.zz = 0.0;
		res.xy = 0.0;
		res.yz = 0.0;
		res.zx = 0.0;
		res.yx = 0.0;
		res.zy = 0.0;
		res.xz = 0.0;
	}

	delete [] d;
	delete [] n;

	return res;
}

///////////////////////////////////////////////////////////////////////////////////////////////
//Hyperfine interactions
void calc_hf_1(int idx, double * hf)
{
    double bx1, by1, bz1;
    double int_const_ne = INT_CONST_NE * RELATIVE_GAMMA_N[idx];

	int i;
	double	x, y, z;

	x = COORDINATES[3 * idx + 0]; y = COORDINATES[3 * idx + 1]; z = COORDINATES[3 * idx + 2];

	double s1x, s1y, s1z, s1n;
	double len3;

	double * d = new double [3];
	double * n = new double [3];
	double len, nx, ny, nz;

	len = 0.;
	for (i = 0; i < 3; i++)
	{
		d[i] = COORDINATES[3 * idx + i] - Qubit_4H_SiC[i];
		len += d[i]*d[i];
	}
	len = sqrt(len);
	for (i = 0; i < 3; i++)	{	n[i] = d[i] / len; }

	s1x = vec_dot(3, S1VALUE, bframeX);	s1y = vec_dot(3, S1VALUE, bframeY);	s1z =vec_dot(3, S1VALUE, bframeZ); s1n =vec_dot(3, S1VALUE, n);
	nx	= vec_dot(3, n, bframeX);		ny = vec_dot(3, n, bframeY);		nz = vec_dot(3, n, bframeZ);


	len3 = len * len * len;
	hf[3*idx + 0] += int_const_ne * ( s1x - 3.0 * s1n * nx ) / len3;
	hf[3*idx + 1] += int_const_ne * ( s1y - 3.0 * s1n * ny ) / len3;
	hf[3*idx + 2] += int_const_ne * ( s1z - 3.0 * s1n * nz ) / len3;

	delete [] d;
	delete [] n;

    int nearQ;
    double magnitude1_old, magnitude1_new;

    if(IsCore)
    {
    	nearQ = isNear(idx);



    	if ( nearQ != -1)
    	{
    		count_core ++;

    		cout << "core spins:" << count_core << " [ " << COORDINATES[3*idx + 0] << ", " << COORDINATES[3*idx + 1] << ", " << COORDINATES[3*idx + 2] << " ]. " << endl;
    		cout << "\t ==>" << "[ " << RELAX_COORD[3*nearQ + 0] << ", " << RELAX_COORD[3*nearQ + 1] << ", " << RELAX_COORD[3*nearQ + 2] << "] " << endl;
    		cout << "diff = " << sqrt(
    				( COORDINATES[3*idx + 0]-RELAX_COORD[3*nearQ + 0])*(COORDINATES[3*idx + 0]-RELAX_COORD[3*nearQ + 0])
    				+(COORDINATES[3*idx + 1]-RELAX_COORD[3*nearQ + 1])*(COORDINATES[3*idx + 1]-RELAX_COORD[3*nearQ + 1])
    				+(COORDINATES[3*idx + 2]-RELAX_COORD[3*nearQ + 2])*(COORDINATES[3*idx + 2]-RELAX_COORD[3*nearQ + 2])) << endl;
    		cout << endl;


    		cout << "\t    " << "[ "<< hf[3*idx + 0] << ", " << hf[3*idx + 1] << ", " << hf[3*idx + 2] << "], " << endl;
    		magnitude1_old = sqrt( hf[3*idx + 0]* hf[3*idx + 0]+ hf[3*idx + 1]* hf[3*idx + 1]+ hf[3*idx + 2]* hf[3*idx + 2]);


    		COORDINATES[3*idx + 0] = RELAX_COORD[3*nearQ + 0];
    		COORDINATES[3*idx + 1] = RELAX_COORD[3*nearQ + 1];
    		COORDINATES[3*idx + 2] = RELAX_COORD[3*nearQ + 2];

    		bx1 = S1VALUE[0]*HF_MAT_FIRSTPRINCIPLE[nearQ][0*3+0] + S1VALUE[1]*HF_MAT_FIRSTPRINCIPLE[nearQ][1*3+0] + S1VALUE[2]*HF_MAT_FIRSTPRINCIPLE[nearQ][2*3+0];
    		by1 = S1VALUE[0]*HF_MAT_FIRSTPRINCIPLE[nearQ][0*3+1] + S1VALUE[1]*HF_MAT_FIRSTPRINCIPLE[nearQ][1*3+1] + S1VALUE[2]*HF_MAT_FIRSTPRINCIPLE[nearQ][2*3+1];
    		bz1 = S1VALUE[0]*HF_MAT_FIRSTPRINCIPLE[nearQ][0*3+2] + S1VALUE[1]*HF_MAT_FIRSTPRINCIPLE[nearQ][1*3+2] + S1VALUE[2]*HF_MAT_FIRSTPRINCIPLE[nearQ][2*3+2];

    		hf[3*idx + 0] = (bx1*bframeX[0] + by1*bframeX[1] + bz1*bframeX[2])*1E03 *2.0*PI * 0.5;
    		hf[3*idx + 1] = (bx1*bframeY[0] + by1*bframeY[1] + bz1*bframeY[2])*1E03 *2.0*PI * 0.5;
    		hf[3*idx + 2] = (bx1*bframeZ[0] + by1*bframeZ[1] + bz1*bframeZ[2])*1E03 *2.0*PI * 0.5;

    		cout << "\t ==>" << "[ "<< hf[3*idx + 0] << ", " << hf[3*idx + 1] << ", " << hf[3*idx + 2] << "], " << endl;

    		magnitude1_new = sqrt( hf[3*idx + 0]* hf[3*idx + 0]+ hf[3*idx + 1]* hf[3*idx + 1]+ hf[3*idx + 2]* hf[3*idx + 2]);

    		cout << endl;
    		cout << "hf1: " << magnitude1_old <<" ==> " << magnitude1_new << endl;
    		cout << endl;
    	}
    }
}

void calc_hf_2(int idx, double * hf)
{
    double bx2, by2, bz2;
    double int_const_ne = INT_CONST_NE * RELATIVE_GAMMA_N[idx];

	int i;
	double	x, y, z;

	x = COORDINATES[3 * idx + 0]; y = COORDINATES[3 * idx + 1]; z = COORDINATES[3 * idx + 2];

	double s2x, s2y, s2z, s2n;
	double len3;

	double * d = new double [3];
	double * n = new double [3];
	double len, nx, ny, nz;

	len = 0.;
	for (i = 0; i < 3; i++)
	{
		d[i] = COORDINATES[3 * idx + i] - Qubit_4H_SiC[i];
		len += d[i]*d[i];
	}
	len = sqrt(len);
	for (i = 0; i < 3; i++)	{	n[i] = d[i] / len; }

	s2x = vec_dot(3, S2VALUE, bframeX);	s2y = vec_dot(3, S2VALUE, bframeY);	s2z =vec_dot(3, S2VALUE, bframeZ); s2n =vec_dot(3, S2VALUE, n);
	nx	= vec_dot(3, n, bframeX);		ny = vec_dot(3, n, bframeY);		nz = vec_dot(3, n, bframeZ);

	len3 = len * len * len;
	hf[3*idx + 0] += int_const_ne * ( s2x - 3.0 * s2n * nx ) / len3;
	hf[3*idx + 1] += int_const_ne * ( s2y - 3.0 * s2n * ny ) / len3;
	hf[3*idx + 2] += int_const_ne * ( s2z - 3.0 * s2n * nz ) / len3;

	delete [] d;
	delete [] n;

    int nearQ;
    double  magnitude2_old,magnitude2_new;
    if(IsCore)
    {
    	nearQ = isNear(idx);
    	if ( nearQ != -1 )
    	{
    		count_core ++;

    		cout << "core spins:" << count_core << " [ " << COORDINATES[3*idx + 0] << ", " << COORDINATES[3*idx + 1] << ", " << COORDINATES[3*idx + 2] << " ]. " << endl;
    		cout << "\t ==>" << "[ " << RELAX_COORD[3*nearQ + 0] << ", " << RELAX_COORD[3*nearQ + 1] << ", " << RELAX_COORD[3*nearQ + 2] << "] " << endl;
    		cout << "diff = " << sqrt(
    				( COORDINATES[3*idx + 0]-RELAX_COORD[3*nearQ + 0])*(COORDINATES[3*idx + 0]-RELAX_COORD[3*nearQ + 0])
    				+(COORDINATES[3*idx + 1]-RELAX_COORD[3*nearQ + 1])*(COORDINATES[3*idx + 1]-RELAX_COORD[3*nearQ + 1])
    				+(COORDINATES[3*idx + 2]-RELAX_COORD[3*nearQ + 2])*(COORDINATES[3*idx + 2]-RELAX_COORD[3*nearQ + 2])) << endl;
    		cout << endl;


    		cout << "\t    " << "[ "<< hf[3*idx + 0] << ", " << hf[3*idx + 1] << ", " << hf[3*idx + 2] << "], " << endl;
    		magnitude2_old = sqrt( hf[3*idx + 0]* hf[3*idx + 0]+ hf[3*idx + 1]* hf[3*idx + 1]+ hf[3*idx + 2]* hf[3*idx + 2]);


    		COORDINATES[3*idx + 0] = RELAX_COORD[3*nearQ + 0];
    		COORDINATES[3*idx + 1] = RELAX_COORD[3*nearQ + 1];
    		COORDINATES[3*idx + 2] = RELAX_COORD[3*nearQ + 2];

    		bx2 = S2VALUE[0]*HF_MAT_FIRSTPRINCIPLE[nearQ][0*3+0] + S2VALUE[1]*HF_MAT_FIRSTPRINCIPLE[nearQ][1*3+0] + S2VALUE[2]*HF_MAT_FIRSTPRINCIPLE[nearQ][2*3+0];
    		by2 = S2VALUE[0]*HF_MAT_FIRSTPRINCIPLE[nearQ][0*3+1] + S2VALUE[1]*HF_MAT_FIRSTPRINCIPLE[nearQ][1*3+1] + S2VALUE[2]*HF_MAT_FIRSTPRINCIPLE[nearQ][2*3+1];
    		bz2 = S2VALUE[0]*HF_MAT_FIRSTPRINCIPLE[nearQ][0*3+2] + S2VALUE[1]*HF_MAT_FIRSTPRINCIPLE[nearQ][1*3+2] + S2VALUE[2]*HF_MAT_FIRSTPRINCIPLE[nearQ][2*3+2];

    		hf[3*idx + 0] = (bx2*bframeX[0] + by2*bframeX[1] + bz2*bframeX[2])*1E03 * 2.0*PI * 0.5;     // the factor 0.5 arises from spin 1/2 of C-13
    		hf[3*idx + 1] = (bx2*bframeY[0] + by2*bframeY[1] + bz2*bframeY[2])*1E03 * 2.0*PI * 0.5;     // convert from MHz to k rad / s
    		hf[3*idx + 2] = (bx2*bframeZ[0] + by2*bframeZ[1] + bz2*bframeZ[2])*1E03 * 2.0*PI * 0.5;

    		cout << "\t ==>" << "[ "<< hf[3*idx + 0] << ", " << hf[3*idx + 1] << ", " << hf[3*idx + 2] << "], " << endl;

    		magnitude2_new = sqrt( hf[3*idx + 0]* hf[3*idx + 0]+ hf[3*idx + 1]* hf[3*idx + 1]+ hf[3*idx + 2]* hf[3*idx + 2]);

    		cout << endl;
    		cout << "; hf2: " << magnitude2_old << " ==> " << magnitude2_new << endl;
    		cout << endl;
    	}
    }
}

///////////////////////////////////////////////////////////////////////////////////////////////
//Hyperfine mediated interactions
void calc_hf_tensor(int idx)
{
    double int_const_ne = INT_CONST_NE * RELATIVE_GAMMA_N[idx];

	int i;
	double	x, y, z;

	x = COORDINATES[3 * idx + 0]; y = COORDINATES[3 * idx + 1]; z = COORDINATES[3 * idx + 2];

	double * d = new double [3];
	double * n = new double [3];
	double len, nx, ny, nz, len3;

	len = 0.;
	for (i = 0; i < 3; i++)
	{
		d[i] = COORDINATES[3 * idx + i] - Qubit_4H_SiC[i];
		len += d[i]*d[i];
	}
	len = sqrt(len);
	len3 = len * len * len;
	for (i = 0; i < 3; i++)	{	n[i] = d[i] / len; }


	nx	= vec_dot(3, n, bframeX);		ny = vec_dot(3, n, bframeY);		nz = vec_dot(3, n, bframeZ);

	HF_TENSOR[idx * 9 + 0] = 1.0 - 3.0 * nx * nx;	HF_TENSOR[idx * 9 + 1] = 	- 3.0 * nx * ny;	HF_TENSOR[idx * 9 + 2] = 	 - 3.0 * nx * nz;
	HF_TENSOR[idx * 9 + 3] = 	 - 3.0 * ny * nx;	HF_TENSOR[idx * 9 + 4] = 1.0- 3.0 * ny * ny;	HF_TENSOR[idx * 9 + 5] = 	 - 3.0 * ny * nz;
	HF_TENSOR[idx * 9 + 6] =	 - 3.0 * nz * nx;	HF_TENSOR[idx * 9 + 7] =    - 3.0 * nz * ny;	HF_TENSOR[idx * 9 + 8] = 1.0 - 3.0 * nz * nz;

	for(i = 0; i < 9; i ++)
		HF_TENSOR[idx * 9 + i] *= (int_const_ne / len3);
}

void calc_mediated_int(int index_i, int index_j, int state, double * interaction)
{
	int i;

	int perturb_state1, perturb_state2;
	if(state == 0) {perturb_state1 = 1; perturb_state2 = 2;}
	if(state == 1) {perturb_state1 = 0; perturb_state2 = 2;}
	if(state == 2) {perturb_state1 = 0; perturb_state2 = 1;}

	int alpha_beta_1, alpha_beta_2, beta_alpha_1, beta_alpha_2;
	alpha_beta_1 = state * 3 + perturb_state1;
	alpha_beta_2 = state * 3 + perturb_state2;
	beta_alpha_1 = perturb_state1 * 3 + state;
	beta_alpha_2 = perturb_state2 * 3 + state;

	double delta_1, delta_2;
	double degenerate_eps = 100.0;// k rad / s
	delta_1 = EIGEN_VALUES[state] - EIGEN_VALUES[perturb_state1];
	delta_2 = EIGEN_VALUES[state] - EIGEN_VALUES[perturb_state2];
	if(fabs(delta_1) < degenerate_eps || fabs(delta_2) < degenerate_eps)
	{
		cout << "delta=" << delta_1 << ", " << delta_2 <<endl;
		assert(fabs(delta_1) > degenerate_eps && fabs(delta_2) > degenerate_eps);
	}

	double * hftensor_i, * hftensor_j, * Rmat1, * Rmat2;
	CLD * S_ab1, * S_ba1, * S_ab2, * S_ba2;
	hftensor_i = HF_TENSOR + index_i * 9;
	hftensor_j = HF_TENSOR + index_j * 9;

	Rmat1 = RMatrix[alpha_beta_1];
	Rmat2 = RMatrix[alpha_beta_2];
	S_ab1 = S_ab_Bframe[alpha_beta_1];
	S_ba1 = S_ab_Bframe[beta_alpha_1];
	S_ab2 = S_ab_Bframe[alpha_beta_2];
	S_ba2 = S_ab_Bframe[beta_alpha_2];

	double * temp = new double [3 * 3];
	double * temp1 = new double [3 * 3];
	double * temp2 = new double [3 * 3];

	char * noTrans;
	noTrans	= "NoTranspose";

	MKL_real_mat_mult(3, noTrans, noTrans, 1.0, Rmat1, hftensor_j, 0.0, temp);
	MKL_real_mat_mult(3, noTrans, noTrans, 1.0, hftensor_i, temp, 0.0, temp1);

	MKL_real_mat_mult(3, noTrans, noTrans, 1.0, Rmat2, hftensor_j, 0.0, temp);
	MKL_real_mat_mult(3, noTrans, noTrans, 1.0, hftensor_i, temp, 0.0, temp2);


	CLD * sab_A_vec_1 = new CLD [3];
	CLD * sba_A_vec_1 = new CLD [3];
	CLD * sab_A_vec_2 = new CLD [3];
	CLD * sba_A_vec_2 = new CLD [3];
	double hf_mediate_x, hf_mediate_y, hf_mediate_z;

	if(index_i != index_j)
	{
		for(i = 0; i < 9; i++)
		{
			if(fabs(delta_1) > degenerate_eps)
				interaction[i] = temp1[i] / delta_1;

			if(fabs(delta_2) > degenerate_eps)
				interaction[i] += temp2[i] / delta_2;
		}
	}
	else
	{
		for(i = 0; i < 9; i++)
			interaction[i] = 0.0;

		vec_dot_mat(S_ab1, hftensor_i, sab_A_vec_1);
		vec_dot_mat(S_ba1, hftensor_i, sba_A_vec_1);
		vec_dot_mat(S_ab2, hftensor_i, sab_A_vec_2);
		vec_dot_mat(S_ba2, hftensor_i, sba_A_vec_2);

		hf_mediate_x = - ( (sab_A_vec_1[1] * sba_A_vec_1[2] / delta_1).imag() + (sab_A_vec_2[1] * sba_A_vec_2[2] / delta_2).imag() );
		hf_mediate_y = - ( (sab_A_vec_1[2] * sba_A_vec_1[0] / delta_1).imag() + (sab_A_vec_2[2] * sba_A_vec_2[0] / delta_2).imag() );
		hf_mediate_z = - ( (sab_A_vec_1[0] * sba_A_vec_1[1] / delta_1).imag() + (sab_A_vec_2[0] * sba_A_vec_2[1] / delta_2).imag() );

		if(state == STATE_1)
		{
			NE_INT_MAT1[3*index_i + 0] += hf_mediate_x;
			NE_INT_MAT1[3*index_i + 1] += hf_mediate_y;
			NE_INT_MAT1[3*index_i + 2] += hf_mediate_z;
		}
		else
		{
			NE_INT_MAT2[3*index_i + 0] += hf_mediate_x;
			NE_INT_MAT2[3*index_i + 1] += hf_mediate_y;
			NE_INT_MAT2[3*index_i + 2] += hf_mediate_z;
		}
	}

	delete [] temp;
	delete [] temp1;
	delete [] temp2;
	delete [] sab_A_vec_1;
	delete [] sba_A_vec_1;
	delete [] sab_A_vec_2;
	delete [] sba_A_vec_2;
}

///////////////////////////////////////////////////////////////////////////////////////////////
//Assistant functions: coordinate/basis transformation
void set_b_frame()
{
	double * z;
	z = new double [3];


	z[0] =  Zdir_4H_SiC[0];
	z[1] =  Zdir_4H_SiC[1];
	z[2] =  Zdir_4H_SiC[2];

	if (THETA == 0.)
	{
		bframeX[0] = Xdir_4H_SiC[0];	bframeX[1] = Xdir_4H_SiC[1];	bframeX[2] = Xdir_4H_SiC[2];
		bframeY[0] = Ydir_4H_SiC[0];	bframeY[1] = Ydir_4H_SiC[1];	bframeY[2] = Ydir_4H_SiC[2];
		bframeZ[0] = Zdir_4H_SiC[0];	bframeZ[1] = Zdir_4H_SiC[1];	bframeZ[2] = Zdir_4H_SiC[2];
	}
	else
	{
		bframeZ[0] = sin(THETA)*cos(PHI)*Xdir_4H_SiC[0] + sin(THETA)*sin(PHI)*Ydir_4H_SiC[0] + cos(THETA) * Zdir_4H_SiC[0];
		bframeZ[1] = sin(THETA)*cos(PHI)*Xdir_4H_SiC[1] + sin(THETA)*sin(PHI)*Ydir_4H_SiC[1] + cos(THETA) * Zdir_4H_SiC[1];
		bframeZ[2] = sin(THETA)*cos(PHI)*Xdir_4H_SiC[2] + sin(THETA)*sin(PHI)*Ydir_4H_SiC[2] + cos(THETA) * Zdir_4H_SiC[2];

		vec3_cross(z, bframeZ, bframeY, true);			// bframeY = cross(z, B)
		vec3_cross(bframeY, bframeZ, bframeX, true);	// bframeX = cross(Y, B)
	}
	delete [] z;

}

void Crystal2Nature(CLD * coord_c, CLD * coord_n)
{
	//coord_n(x,y,z) = coord_c(x,y,z) . matrix(row_1: Xdir_4H_SiC, row_2: Ydir_4H_SiC, row_3: Zdir_4H_SiC)
	coord_n[0] = coord_c[0] * Xdir_4H_SiC[0] + coord_c[1] * Ydir_4H_SiC[0] + coord_c[2] * Zdir_4H_SiC[0];
	coord_n[1] = coord_c[0] * Xdir_4H_SiC[1] + coord_c[1] * Ydir_4H_SiC[1] + coord_c[2] * Zdir_4H_SiC[1];
	coord_n[2] = coord_c[0] * Xdir_4H_SiC[2] + coord_c[1] * Ydir_4H_SiC[2] + coord_c[2] * Zdir_4H_SiC[2];
}

void Nature2Bframe(CLD * coord_n, CLD * coord_b)
{
	//coord_b(x,y,z) = ( coord_c(x,y,z) . bframeX, coord_c(x,y,z) . bframeY, coord_c(x,y,z) . bframeZ )
	coord_b[0] = coord_n[0] * bframeX[0] + coord_n[1] * bframeX[1] + coord_n[2] * bframeX[2];
	coord_b[1] = coord_n[0] * bframeY[0] + coord_n[1] * bframeY[1] + coord_n[2] * bframeY[2];
	coord_b[2] = coord_n[0] * bframeZ[0] + coord_n[1] * bframeZ[1] + coord_n[2] * bframeZ[2];
}

void uni_transformation(int dim, CLD * mat, CLD * R_mat, CLD * res_mat)
{
	CLD * temp = new CLD [dim * dim];

    char * noTrans, * conjTrans, * Trans;

    conjTrans	= "ConjugateTranspose";
    noTrans		= "NoTranspose";
    Trans		= "Transpose";		//Transpose Spin matrix Sx, Sy, and Sz according to the MKL matrix storage convention.

	MKL_mat_mult(3, Trans, noTrans, COMPLEX_ONE, mat, R_mat, COMPLEX_ZERO, temp);

	MKL_mat_mult(3, conjTrans, noTrans, COMPLEX_ONE, R_mat, temp, COMPLEX_ZERO, res_mat);

	delete [] temp;
}

///////////////////////////////////////////////////////////////////////////////////////////////
//Assistant functions: simple linear algebar
double vec_dot(int dim, double * vec1, double * vec2)
{
    double res;
    int i;

    res = 0.0;
    for (i = 0; i < dim; i++)
    {
            res += vec1[i]*vec2[i];
    }

    return(res);
}

void vec3_cross(double * vec1, double * vec2, double * vec3, bool is_norm)
{
    double len;
    vec3[0] = vec1[1]*vec2[2] - vec1[2]*vec2[1];
    vec3[1] = vec1[2]*vec2[0] - vec1[0]*vec2[2];
    vec3[2] = vec1[0]*vec2[1] - vec1[1]*vec2[0];
    if(is_norm)
    {
            len = sqrt(vec3[0]*vec3[0] + vec3[1]*vec3[1] + vec3[2]*vec3[2]);
            vec3[0] /= len; vec3[1] /= len; vec3[2] /= len;
    }
}

void vec_dot_mat(CLD * S_ab, double * hftensor_i, CLD * sab_A_vec)
{
	sab_A_vec[0] = S_ab[0] * hftensor_i[0] + S_ab[1] * hftensor_i[1] +S_ab[2] * hftensor_i[2];
	sab_A_vec[1] = S_ab[0] * hftensor_i[3] + S_ab[1] * hftensor_i[4] +S_ab[2] * hftensor_i[5];
	sab_A_vec[2] = S_ab[0] * hftensor_i[6] + S_ab[1] * hftensor_i[7] +S_ab[2] * hftensor_i[8];
}
